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Abstract. In this paper, we derive a variational integrator for certain higlily oscillatory problems 
in mechanics. To do this, we take a new approach to the splitting of fast and slow potential forces: 
rather than splitting these forces at the level of the differential equations or the Hamiltonian, 
we split the two potentials with respect to the Lagrangian action integral. By using a different 
quadrature rule to approximate the contribution of each potential to the action, we arrive at a 
geometric integrator that is implicit in the fast force and explicit in the slow force. This can allow 
for significantly longer time steps to be taken (compared to standard explicit methods, such as 
Stormer/Verlet) at the cost of only a linear solve rather than a full nonlinear solve. We also analyze 
the stability of this method, in particular proving that it eliminates the linear resonance instabilities 
that can arise with explicit multiple-time-stepping methods. Next, we perform some numerical 
experiments, studying the behavior of this integrator for two test problems: a system of coupled 
linear oscillators, for which we compare against the resonance behavior of the r-RESPA method; 
and slow energy exchange in the Fermi-Pasta-Ulam problem, which couples fast linear oscillators 
with slow nonlinear oscillators. Finally, we prove that this integrator accurately preserves the slow 
energy exchange between the fast oscillatory components, which explains the numerical behavior 
observed for the Fermi-Pasta-Ulam problem. 



1. Introduction 



1.1. Problem Background. Many systems in Lagrangian mechanics have components acting on 
different time scales, posing a challenge for traditional numerical integrators. Examples include: 

(1) Elasticity: Several spatial elements of varying stiffness, resulting from irregular meshes 
and/or inhomogeneous materials |Lewetal^|2003| . 

(2) Planetary Dynamics: A/^-body problem with nonlinear gravitational forces, arising from 
pairwise inverse-square potentials. Multiple time scales result from the different distances 



between the bodies i Fair and Bertschinger 2007 . 
(3) Highly Oscillatory Problems: Potential energy can be split into a "fast" linear oscillatory 
component and a "slow" nonlinear component. These problems are widely encountered 
in modeling molecular dynamics jLeimkuhler et al.j 1996| , but have also been used to 



model other diverse applications, for example, in computer animation fEberhardt et al.| 
[20001 [Boxerman and Ascherl [2004t . 
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Because these systems each satisfy a Lagrangian variational principle, they lend themselves readily 
to variational integrators: a class of geometric numerical integrators designed for simulating 
Lagrangian mechanical systems. By construction, variational integrators preserve a discrete 
version of this Lagrangian variational structure; consequently, they are automatically symplectic 



and momentum- conserving, with good long-time energy behavior ^Marsden and West 2001 



Explicit Methods, Multiple Time Stepping, and Resonance Instability. The Stormer/Verlet (or 
leapfrog) method is one of the canonical examples of a geometric (and variational) numerical 



integrator (see Hairer et al. 2003^ . Yet, it and other simple, explicit time stepping methods do not 



perform well for problems with multiple time scales. The maximum stable time step for these 
methods is dictated by the stiffest mode of the underlying system; therefore, the fastest force 
dictates the number of evaluations that must be taken for all forces, despite the fact that the 
slow-scale forces may be (and often are) much more expensive to evaluate. 

To reduce the number of costly function evaluations associated to the slow force, several explicit 
variational integrators use multiple time stepping, whereby different time step sizes are used to 
advance the fast and slow degrees of freedom. These include substepping methods, such as 
Verlet-I/r-RESPA and mollified impulse, where for each slow time step, an integer number of fast 
substeps are taken jlzaguirre et al.[ 2002 . More recently, asynchronous variational integrators 



(AVIs) have been developed, removing the restriction for fast and slow time steps to be integer (or 
even rational) multiples of one another jLew et aL| 2003| . Multiple-time-stepping methods can 



be more efficient than single-time-stepping explicit methods, like Stormer/Verlet, since one can 
fully resolve the fast oscillations while taking many fewer evaluations of the slow forces. This is 
especially advantageous for highly oscillatory problems, where the slow forces are nonlinear and 
hence more computationally expensive to evaluate. 

One drawback of multiple-time-stepping methods, however, is that they can exhibit linear 
resonance instability. This phenomenon occurs when the slow impulses are nearly synchronized, 
in phase, with the the fast oscillations. These impulses artificially drive the system at a resonant 
frequency, causing the energy (and hence the numerical error) to increase without bound. The 



problem of numerical resonance is well knovm for substepping methods Biesiadecki and Skeel 



1993 1, and has also recently been shown for AVIs as well — in fact, the subset of fast and slow time 



step size pairs leading to resonance instability is dense in the space of all possible parameters i Fong 



et al. 2007 i. Resonance instability can therefore be difficult to avoid, particularly in highly 



oscillatory systems with many degrees of freedom, as in molecular dynamics applications. 

Implicit Methods for Single Time Stepping with Longer Step Sizes. Because multiple-time-stepping 
methods have these resonance problems, a number of single-time-stepping methods have been 
developed specifically for highly oscillatory problems. As noted earlier, single-time-stepping 
methods cannot fully resolve the fast oscillations without serious losses in efficiency. Therefore, 
the goal of these methods is to take long time steps, without actually resolving the fast oscillations, 
while still accurately capturing the macroscopic behavior that emerges from the coupling between 
fast and slow scales. The challenge is to design methods that allow for these longer time steps, 
without destroying either numerical stability or geometric structure. 
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One obvious candidate integrator is tiie implicit midpoint method, wiiich is (linearly) uncon- 
ditionally stable, as well as variational (hence symplectic) and symmetric. Unfortunately, the 
stability of the method comes at a cost: because the integrator is implicit in the slow force, which is 
generally nonlinear, a nonlinear system of equations must be solved at every time step. Therefore, 
just like the fully-resolved Stormer/Verlet method, this means that the implicit midpoint method 
requires an excessive number of function evaluations. 

Implicit- Explicit Integration. For highly oscillatory problems, implicit-explicit (IMEX) integrators 
have been proposed as a potentially attractive alternative to either explicit, multiple-time-stepping 
methods or implicit, single-time-stepping methods. Rather than using separate fast and slow time 
step sizes, IMEX methods combine implicit integration (e.g., backward Euler) for the fast force 
with explicit integration (e.g., forward Euler) for the slow force. Because the fast force is linear, 
this semi-implicit approach requires only a linear solve for the implicit portion, as opposed to the 
expensive nonlinear solve that would be required for a fully implicit integrator, like the implicit 
midpoint method. 



IMEX methods were developed by Crouzeix 1980 , and have continued to progress, including 



the introduction of IMEX Runge-Kutta schemes for PDEs by Ascher et al. 1997 . However, in all of 
these methods, the splitting is done at the level of the Euler-Lagrange differential equations, rather 
than at the variational level of the Lagrangian. Consequently, a wide variety of IMEX schemes 
have been created, both geometric and non-geometric, but in general they cannot guarantee 
properties such as symplecticity, momentum conservation, or good long-time energy behavior, 
which automatically hold for variational integrators. As an example of an IMEX integrator that 



is not "geometric" in the usual sense, consider the LI and LIN methods of Zhang and Schlick 



1993 1, which combine the backward Euler method with explicit Langevin dynamics for molecular 
systems. In particular, to ameliorate the artificial numerical dissipation introduced by using 
backward Euler, these methods rely on stochastic forcing to inject the missing energy back into 
the system. 

In this paper, we develop IMEX numerical integration from a Lagrangian, variational point 
of view. We do this by splitting the fast and slow potentials at the level of the Lagrangian action 
integral, rather than with respect to the differential equations or the Hamiltonian. From this 
viewpoint, implicit- explicit integration is an automatic consequence of discretizing the action 
integral using two distinct quadrature rules for the slow and fast potentials. The resulting discrete 
Euler-Lagrange equations coincide with a semi-implicit algorithm that was originally introduced 
by jZhang and Skeel 1997| as a "cheaper" alternative to the implicit midpoint method; Ascher 



and Reich] jl999b I also studied a variant of this method for certain problems in molecular dy- 
namics, replacing the implicit midpoint step by the energy- conserving (but non-symplectic) 
Simo-Gonzales method. 

We also show that this variational IMEX method is free of resonance instabilities; the proof of 
this fact is naturally developed at the level of the Lagrangian, and does not require an examination 
of the associated Euler-Lagrange equations. We then compare the resonance-free behavior of 
variational IMEX to the multiple-time-stepping method r-RESPA in a numerical simulation of 
coupled slow and fast oscillators. Next, we evaluate the stability of the variational IMEX method. 
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for large time steps, in a computation of slow energy exchange in the Fermi-Pasta-Ulam problem. 
Finally, we prove that the variational IMEX method accurately preserves this slow energy exchange 
behavior (as observed in the numerical experiments) by showing that it corresponds to a modified 
impulse method. 



1.2. A Brief Review of Variational Integrators. The idea of variational integrators was studied 
by Suris |1990| and Moser and Veselov'fl991 , among others, and a general theory was developed 



over the subsequent decade (see Marsden and West||2001l for a comprehensive survey). 



Suppose we have a mechanical system on a configuration manifold Q, specified by a Lagrangian 
L: TQ M. Given a set of discrete time points to < •■• < with uniform step size h, we wish 
to compute a numerical approximation q„ qitn), n = 0,...,N, to the continuous trajectory 
cj{t). To construct a variational integrator for this problem, we define a discrete Lagrangian 
Lh'. Q X Q ^M., replacing tangent vectors by pairs of consecutive configuration points, so that 
with respect to some interpolation method and numerical quadrature rule we have 



H + l 



Lh{qn,qn+i]^ \ L[q,q)dt. 



Then the action integral over the whole time interval is approximated by the discrete action sum 



A/-1 ntN 
n=0 J to 



'to 

If we apply Hamilton's principle to this action sum, so that 5Sh[q] = when variations are 
taken over paths with fixed endpoints, then this yields the discrete Euler-Lagrange equations 

DiLh {qn,qn+i) +£'2^/2 (^7«-i,^n) =0, n = \,...,N -\, 

where Di and D2 denote partial differentiation in the first and second arguments, respectively. 
This defines a two-step numerical method on Q x Q, mapping [qn-iy^n] * (^7n>^n+i)- The 
equivalent one-step method on the cotangent bundle T*Q, mapping [qn,Pn) ^ . > is 
defined by the discrete Legendre transform 

Pn = -DiLh [qnyCln+l] > Pn+1 = D2Lh (<7«,<7«+l) , 

where the first equation updates q, and the second updates p. 



Examples. Consider a Lagrangian of the form L [q, q] — ^q^Mq — V[q), where Q = M'', M is a 
constant d x d mass matrix, and 1^: Q — > M is a potential. If we use linear interpolation of q with 
trapezoidal quadrature to approximate the contribution of V to the action integral, we get 

Ttrap f \ _ ^ f ^ri+l ~ qn\^ f qn+1 ~ qn \ , Vicin) + V[q„+i) 

[q„,q„+,]--]^ - J M]^ - j-h - , 

which we call the trapezoidal discrete Lagrangian. It is straightforward to see that the discrete 
Euler-Lagrange equations for l}^^ correspond to the explicit Stormer/Verlet method. Alterna- 
tively, if we use midpoint quadrature to approximate the integral of the potential, this yields the 
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midpoint discrete Lagrangian, 

for which the resulting integrator is the implicit midpoint method. 



2. A Variational IMEX Method 

In this section, we show how to develop a variational integrator that combines aspects of the 
Stormer/Verlet and implicit midpoint methods mentioned above. The main idea is that, given a 
splitting of the potential energy into fast and slow components, we define the discrete Lagrangian 
by applying the midpoint quadrature rule to the fast potential and the trapezoidal quadrature rule 
to the slow potential. The resulting variational integrator is implicit in the fast force and explicit 
in the slow force. After this, we focus on the specific case of highly oscillatory problems, where 
the fast potential is quadratic (corresponding to a linear fast force). In this case, we show that the 
IMEX integrator can be understood as Stormer/Verlet with a modified mass matrix. 

It should be noted that the results in this section can also be verified directly in terms of 
the numerical algorithm, and do not strictly require making use of the Lagrangian variational 
structure. However, we find the variational perspective to be useful and illustrative, both in 
arriving at this particular IMEX algorithm and in interpreting its numerical features. 



2.1. The IMEX Discrete Lagrangian and Equations of Motion. Suppose that we have a Lagran- 
gian of the form L [q, q) = ^q^Mq — U{q) — W{q), where f/ is a slow potential and 14^ is a fast 
potential, for the configuration space Q = . Then define the IMEX discrete Lagrangian 



h ( qn+l -(^riY ^n+l - qn\ , U {qn) + U {qn+l) 



M 



+ qn+l j 



-,"-.7 2\ h ) "\ h J ' 2 - \ 2 

using (explicit) trapezoidal approximation for the slow potential and (implicit) midpoint approxi- 
mation for the fast potential. The discrete Euler-Lagrange equations give the two-step variational 
integrator on Q x Q 



qn+i-2q„ + qn-i = -h^M ^ 
and the corresponding discrete Legendre transform is given by 

To see how this translates into an algorithm for a one-step integrator on T*Q, it is helpful to 
introduce the intermediate stages 



+ <7n+l A 
2 )' 



Pt = Pn-\yU{qn), 



P„+l = Pn+l+ 2^^Wn+0' 
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Substituting these into the previous expression and rearranging yields the algorithm 

, h . . 

Stepl: Pn=Pn--'^U[qn), 



qn+i =qn + hM ^ 



Step 2: 



Pt + Pn+l 

2 



+ ,^...fqn + qn+i\ 

Pn+l=Pt-h^^[ 2 J' 



Step 3: Pn+i = p„+i--'^U[qn+i), 

where Step 2 corresponds to a step of the implicit midpoint method. 
This can be summarized, in the style of impulse methods, as: 

(1) kick: explicit kick from f/ advances {qn,Pn) {j^n,pX), 

(2) oscillate: implicit midpoint method with Vl^ advances {qn,pX) ^ (^n+i.P^+i). 

(3) kick: explicit kick from f7 advances {^qn+i,p~n+i} ^ {qn+i,Pn+i]- 

In particular, notice that this reduces to the Stormer/Verlet method when VW = and to the 
implicit midpoint method when Vf/ = 0. Also, if the momentum p„ does not actually need to be 
recorded at the full time step (i.e., collocated with the position then Step 3 can be combined 
with Step 1 of the next iteration to create a staggered "leapfrog" method. 



Interpretation as a Hamiltonian Splitting Method. This algorithm on T*Q can also be interpreted 
as a fast-slow splitting method i jMcLachlan and Quispel[|2002[|Hairer et aL||2006| II.5 and VHI.4.1) 
for the separable Hamiltonian H = T +U +W, where T is the kinetic energy, as follows. Let 
^r+w. j*Q _^ j,*Q (jgjjQi^g ^^j^g numerical flow of the implicit midpoint method with time step 
size h, applied to the fast portion of the Hamiltonian T+W, and let ip^: T*Q — > T*Q be the exact 
Hamiltonian flow for the slow potential U (i.e., constant acceleration without displacement). 
Then the variational IMEX method has the flow map '■ T*Q — > T*Q, which can be written as the 
following composition of exact and numerical flows: 

This formulation highlights the fact that variational IMEX is symmetric (since it is a symmetric 
composition of symmetric methods) as well as symplectic (since it can be written as a composition 
of symplectic maps). 



2.2. Application to Highly Oscillatory Problems. For highly oscillatory problems on Q = M*^ , we 
start by taking a quadratic fast potential 

W{q) = ^q^rt^q, fl E M^^^ symmetric and positive semidefinite. 

A prototypical n is given by the block-diagonal matrix f2 = ( o ) , where some of the degrees of 
freedom are subjected to an oscillatory force with constant fast frequency £0 » 1. We also denote 
the slow force g{q) = —VU{q) and assume, without loss of generality, that the constant mass 
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matrix is given by M = /. Therefore, the nonlinear system we wish to approximate numerically is 

This is the conventional setup for highly oscillatory problems, used by |Hairer et al. l |2006 XIII) 
and others. 

Applying the IMEX method to this example, we get the discrete Lagrangian 

rlMEXf^ ^ ( ^ri+\ -qn Y ( C]n+1 -Cjn ^ 

_ j^ U{qn) + U{qn+l] __ ^ ^ qn + qn+1 ^ ^2 + 

and so the two-step IMEX scheme is given by the discrete Euler-Lagrange equations 

^7„+l -2^7„ + ^7„_l + — n [qn+i + 2qn + qn-i] = h g(qn)- 
Combining terms, we can rewrite this as 



[qn+i - 2qn + qn-i) + ^ n = ^ g [qn] , 



which is equivalent to Stormer/Verlet with a modified mass matrix / + [h^/lY. This equivalence 
can similarly be shown to hold for the one-step formulation of the IMEX scheme on T*Q — that is, 
the two methods also produce the same pn, as well as the same qn- 

In fact, this correspondence between IMEX and a modified Stormer/Verlet method is true not 
just in the discrete Euler-Lagrange equations, but in the discrete Lagrangian itself. This follows 
immediately from the following proposition. 

Proposition 2.1. Suppose we have a Lagrangian L[q,q) — ^q^Mq — ^q^Q^q and its corre- 
sponding midpoint discrete Lagrangian L™'^. Next, define the modified Lagrangian L [q, q] — 
^q^Mq — ^q^rfiq, having the same quadratic potential but a different mass matrix M , andtake 
its trapezoidal discrete Lagrangian L^!"^. Then L'T^'^ = V["^ when M = M + [hD./2)^. 



Proof. The midpoint discrete Lagrangian is given by 



Now, notice that we can rearrange the terms 



"I 2 J " I 2 J = l 2 J " I 2 )-2'^n^^n--q„^,nqn+. 

Therefore the discrete Lagrangian can be written in the trapezoidal form 

rmidf A h fq„+i-q„Y\,. , f^^Y] fqn+l-qn\ hfl Tr.2 Y T n,2 ^ 



which is precisely V^'^^ {qn,qn+i] when M = M + {hQ./2)'^ 



□ 
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Corollary 2.2. Consider a highly oscillatory system with an arbitrary slow potential U, quadratic 
fast potential W[q) = ^q^rfiq, and constant mass matrix M — I, so that the Lagrangian L and 
IMEX discrete Lagrangian are defined as above. Next, take the modified Lagrangian L with 
the same potentials but different mass matrix M. Then L™^= Zj^™^ when M = I + [hD./2)^. 



2.3. Analysis of Linear Resonance Stability. To study the linear resonance stability of this IMEX 
integrator, we consider a model problem where U and W both correspond to linear oscillators. Let 
U{q) =^q^q and W{q) = ^q^ffiq, where f2 — col for some £0 » 1, and again let the mass matrix 
M = L Although this is something of a "toy problem" — obviously, one could simply combine U 
and W into a single quadratic potential | (1 + co^) q^q — it is illustrative for studying the numerical 
resonance of multiple-time-stepping methods, since the system has no external forcing terms 
and hence no real physical resonance. 

To prove that the IMEX method does not exhibit linear resonance instability, we show that the 
stability condition only requires that the time step be stable for the explicit slow force, and is 



independent of the fast frequency co. The idea of the proof is to use the results from Section 2.2 



showing that the IMEX method is equivalent to Stormer/Verlet with a modified mass matrix, and 
then to apply the well-known stability criteria for Stormer/Verlet. 

In particular, for a harmonic oscillator with unit mass and frequency v, the Stormer/Verlet 
method is linearly stable if and only if |/jv| < 2, as can be shown by a straightforward calculation 
of the eigenvalues of the propagation matrix Hairer et al. 2006[ p. 23). For a system with constant 
mass m and spring constant v^, this condition generalizes to h'^v^ < Am. 



Theorem 2.3. The IMEX method is linearly stable, for the system described above, if and only if 
h<2 (i.e., if and only ifh is a stable time step size for the slow oscillator alone). 



Proof As proved in the previous section, the IMEX method for this system is equivalent to 
Stormer/Verlet with the modified mass matrix / -I- (/zO/2)2. Now, this modified oscillatory system 
has constant mass m — 1 + [hco/lf' and spring constant v'^ — \ + cci^. Therefore, the necessary 
and sufficient condition for linear stability is 

/z2(l + £o2) <4( 1 + — <y2 

and since the h'^co'^ terms cancel on both sides, this is equivalent to h^ <A,oxh<2. □ 



This shows that, in contrast to multiple-time-stepping methods, the IMEX method does not 
exhibit linear resonance instability. In particular, one can interpret the modified mass matrix 
as giving the system an effective frequency of (l + a)^) /(^l + (/jco/2)^j, which attenuates the 
destabilizing high frequencies in the original system. It should be noted that nonZmear instability 
is known to be possible for the implicit midpoint method, although even that can be avoided with 
a time step size restriction that is considerably weaker than that required for explicit methods (see 
lAscher and Reich|[l999al . 
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Figure 1: Maximum energy error of r-RESPA and variational IMEX, integrated over the time 
interval [0, 1000] for a range of parameters o). The r-RESPA method exhibits resonance instability 
near integer values of coh/n, while the variational IMEX method remains stable. 



3. Numerical Experiments 

3.1. Coupled Linear Oscillators. To illustrate the numerical resonance behavior of the varia- 
tional IMEX scheme, as compared with a multiple-time-stepping method, we consider the model 



problem from Section 2.3 for dimension d = I (i.e., Q = R). Figure 1 shows a log plot of the 
maximum absolute error in total energy (i.e., the Hamiltonian) for both r-RESPA and the varia- 
tional IMEX method, for a range of frequencies o). MATLAB simulations were performed over 
the time interval [0, 1000], with fixed time step size h = 0.l, and with the normalized frequency 
coh/n ranging over (0,4.5]. Additionally, to fully resolve the fast oscillations, r-RESPA took 100 fast 
substeps of size h/100 = 0.001 for each full time step of size h. 

The r-RESPA method exhibits "spikes" in the total energy error near integer values of ojh/n, 
corresponding to the parameters where resonance instability develops and the numerical solution 
becomes unbounded. (The finite size of these spikes is due to the fact that the numerical simula- 
tion was run only for a finite interval of time. Interestingly, one also sees "negative spikes," where 
the fast and slow oscillations are exactly out-of-phase and cancel one another.) It should be noted 
that the small substep size of r-RESPA is sufficient for stable integration of the fast force alone; it 
is only the introduction of the slow force that makes things unstable. By contrast, the maximum 
energy error for the variational IMEX method is nearly constant for all values of oj, showing no 



sign of resonance. This is fully consistent with the theoretical result obtained in Theorem 2.3 



3.2. The Fermi-Pasta-Ulam Problem. As an example of a nontrivial highly oscillatory problem 
with nonlinear slow potential, we chose the modified Fermi-Pasta-Ulam (FPU) problem consid- 
ered by jHairer et al. \2006 1.5 and XIII), whose treatment we will now briefly review. The FPU 
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problem consists of 21 unit point masses, which are chained together, in series, by alternating 



weak nonlinear springs and stiff linear springs. (This particular setup is due to Galgani et al. 



1992 and is a variant of the problem originally introduced by Fermi et al.[[T955 ) Clearly, this 



system becomes rather trivial if we make the stiff springs "infinitely stiff," replacing them by rigid 
constraints (as done by some numerical methods, such as SHAKE/RATTLE). However, for finite 
stiffness, the FPU system exhibits interesting dynamics due to the coupling between fast and slow 
springs. 

Let us denote the displacements of the point masses \yj q\,...,qn^^ (where the endpoints 
^0 = (]2i+\ = are taken to be fixed), and their conjugate momenta by pi = cji for i = l,...,2£.ln 
these variables, the FPU system has the Hamiltonian 

1 ^ ^2 ^ i 

^(^' = 2 S (^2/-l + Pli] + X S " ^2i-if - qiif, 

! = 1 / = 1 /=0 

which contains a quadratic potential for the I stiff linear springs, each with frequency a>, and a 
quartic potential for the £ + 1 soft nonlinear (cubic) springs. However, it is helpful to perform the 
coordinate transformation (following Hairer et al. 2006 p. 22) 

qii + <72i-i q2i — q2i-i 



72 ' 

P2i + P2i-l 



Xl,i 



yi,i 



72 ' 

P2i - P2i-l 



V2 72 
so that (modulo rescaling) Xq,/ corresponds to the location of the zth stiff spring's center, Xi,/ corre- 
sponds to its length, and ya,i,yi,i are the respective conjugate momenta. Writing the Hamiltonian 
in these new variables, we have 



I ^ 2 ^ 



! = 1 



+ ■ 



e-1 



i=l 



which considerably simplifies the form of the fast quadratic potential. 

Following the example treated numerically by]Hairer et al. l |2006| ; McLachlan and O'Neale 



2007 , we consider an instance of the FPU problem, integrated over the time interval [0, 200], 
with parameters £ = 3, o) = 50, whose initial conditions are 

xo.M=i, yo,i(o) = i, xi,i(0) = w-i, yi,i(0)=i, 

with zero for all other initial values. This displays an interesting and complex property of the FPU 
problem, called slow energy exchange, which results from the slow nonlinear coupling between 
the stiff springs. If we consider only the energies in the stiff springs, written as 



then the initial conditions start with all of the energy in /i and none in h.h- Over the course of the 
time interval, this energy is transferred in a characteristic way from h to h, gradually transitioning 
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h + l2+ h 




\J2 / 






(a) Reference solution: (b) Stormer/Verlet with h = 

Stormer/Verlet with time step 0.01. 
size fc = 0.001. 



(c) Stormer/Verlet with h ■ 
0.03. 




(d)IMEXwithh = 0.03. 



(e) IMEXwith/i^O.l. 



(f) IMEXwith/z = 0.15. 




(g) IMEXwith/z = 0.2. 



(h)IMEX with;; = 0.25. 



(i)IMEXwith/! = 0.3. 



Figure 2: The IMEX method robustly captures slow energy exchange in the Fermi-Pasta-Ulam 
problem with o) = 50, even for large time steps. Because the fast force is integrated implicitly, 
IMEX remains stable and degrades gradually as the time step size increases — unlike the fully 
explicit Stormer/Verlet method, which rapidly becomes unstable. 



through the middle spring I2. Furthermore, the total stiff energy 1 — + 12 + H remains nearly 
constant, i.e., is an adiabatic invariant of the system. 



Figure 2 shows several numerical simulations of this FPU energy exchange, computed both 



with Stormer/Verlet and with the variational IMEX method, for different choices of time step 
size. The first plot is a reference solution, computed using Stormer/Verlet with h = 0.001, fully 
resolving the fast oscillations. However, we see that the Stormer/Verlet solution's quality and 
stability degrade rapidly as we increase the step size (for h = 0.03, we have hco = 1.5, which is near 
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(a)IMEXwithfc = 0.1 (b) IMEX with /i = 0.3 

Figure 3: Numerical simulation of the FPU problem for T = 4000, which shows the behavior of the 
IMEX method on the oj^ scale. For h = 0.1, we already have hco = 5, yet the oscillatory behavior 
and adiabatic invariant are qualitatively correct. By contrast, for h = 0.3, the method has begun to 
blow up; oscillatory coupling is a drawback of implicit midpoint methods for large time steps. 



the upper end of the stability region \ h(o\ < 2). By contrast, the variational IMEX method performs 
extremely well for h = 0.03-0.15, degrading gradually as the time step size increases. Even as 
the numerical solution begins to undergo serious degradation for h = 0.2-0.3, the qualitative 



structure of the energy exchange behavior between h,l2, h is still maintained. Compare Hairer 



et al. 2006 p. 24, Figure 5.3; see also McLachlan and O'Neale 2007 who examine a wide variety 



of geometric integrators, particularly trigonometric integrators, for the FPU problem, with respect 
to both resonance stability and slow energy exchange. In particular, these authors found that the 
existing trigonometric integrators exhibit a trade-off between correct energy exchange behavior 
and resonance stability, and that these features tend to be mutually exclusive. 

In Figure 3| we show the numerical behavior of the variational IMEX method, applied to this 
same FPU problem, on a longer time scale (T — 4000) and for large time steps — 0.1,0.3). At 
h — Q.\, the IMEX simulation still displays the correct qualitative energy behavior, with respect to 
both the slow energy exchange and the adiabatic invariant /, and the numerical solution remains 
bounded. However, by ^ = 0.3, numerical stability has broken down, as oscillatory coupling in 
the fast modes leads to unbounded amplitude growth. This illustrates one of the drawbacks of 
implicit midpoint- type methods: despite the lack of linear resonances, numerical instability can 
still result for very large time steps due to nonlinear coupling Ascher and Reich[ 1999a|b| . 

This example was chosen to demonstrate that the variational IMEX method does not attain its 
stability merely by "smoothing out" the fast frequencies, in a way that might destroy the structure 
of any fast-slow nonlinear coupling. Rather, despite the fact that it does not resolve the fast 
frequencies, the method is still capable of capturing the complex multiscale interactions seen in 
the FPU problem. 
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4. Analysis of Slow Energy Exchange in the IMEX Method 

In the previous section, the numerical experiments for the Fermi-Pasta-Ulam problem seemed 
to suggest that the variational IMEX method preserves the slow energy exchange between the fast 
oscillatory modes. This is somewhat surprising, since the method does not actually resolve these 
fast oscillations. However, in this section, we will prove that, in fact, this method does accurately 
reproduce the slow energy exchange behavior, as long as the numerical solutions remain bounded. 
This is demonstrated by showing that the variational IMEX method can be understood as a modi- 
fled impulse method; that is, the midpoint step exactly resolves the oscillations of some modified 
differential equation. We can then apply some of the existing theory about numerical energy 
exchange for impulse methods. (It should be noted that impulse methods, which originated with 
the work of |Deuflhard[ 1979[ can be understood as a special case of trigonometric integrators 



when applied to highly oscillatory problems.) 

First, let us rewrite the fast oscillatory system q + rfiq = as the first-order system 

nq\ _f n\ fnq 
P o)[p 

so it follows that the exact solution satisfies 

'^^lq{t + h)\ f cos{hQ.) sin[hQ)\ ( Qq[t) 



p[t + h) J y-sm[hn) cos{hn)j y p[t) 

We will now show that the implicit midpoint method effectively replaces this rotation matrix for 
n by the rotation matrix corresponding to a modified f2. In the transformed coordinates just 
introduced, the implicit midpoint method has the expression 

I -hn/2\ fnqn+i\ f I hn/2\(nqn 



hn/2 I J y pn+i ) y-hn/2 i j {pn 

Therefore, if we take the skew matrix 



A = 



f2 

-n 



it follows that 

'nq, 



^n+A ^^^2)-i (7+ hA/2) (^'^" 

Pn+l J \Pn 



Notice that the expression (/ - hA/2)~^ [I + hA/2) — cayihA) is the Cayley transform, which maps 
skew matrices to special orthogonal matrices (and can be seen as an approximation to the matrix 
exponential map, which gives the exact solution). Hence the stability matrix is special orthogonal, 
so we can write 

'^nq„+i\ f cos(^/jf2j sin (^hQ^\ fnqn 
Pn+i ) \-sin[hn^ cos(hfi)j^p„ 

for some modified frequency f2. Therefore, the stability matrix for the implicit midpoint method 
corresponds to the exac? flow matrix for a modified oscillatory system. 
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As an example, suppose we have Q— (o for some constant frequency co. Applying the 
Cayley transform, it can be seen that the modified frequency o) satisfies 



hco/2 = tan{ha)/2). 



Squaring both sides, this becomes 



[hco/2f = tan^iha)/2)- 



1 — COS(/Z(W) 



l + cos[hd>) 

which finally gives the solution for the modified frequency, 

1 f l-ihoj/2f 

cb = — arccos ^ 

h yi + {hcL>/2f 

Remark. This perspective provides another explanation as to why the variational IMEX method 
does not exhibit resonance: we always have hco < n. In fact, the Cayley transform does not map 
to a rotation by k, except in the limit as hco — > oo. Therefore, for any finite h and o), we will never 
encounter the resonance points corresponding to integer multiples of n. 

As an aside, this also leads to another possible interpretation for the onset of nonlinear in- 



stability, if the time step size h becomes too large (as we saw in Figure 3 1. Since co < n/h, the 
modified frequency co must shrink as h grows. Informally, then, if o) is very small, this can be 
seen as leading to amplitude growth in the fast modes, since it requires less energy to induce this 
amplification. 



Since the implicit midpoint method has now been seen as the exact solution of a modified 
system, we can write the variational IMEX method as the following modified impulse scheme: 

h 

Pn 



Step 1: 
Step 2: 
Step 3: 



= Pn--yU[q„), 



P'n+l 



COS 

-sin 



sin 
cos 



P«+i = P„+i 



-VU{qn+i). 



VLq„ 
n+ 

r n 



Suppose again that Q— (o for some constant frequency co, and likewise ^=[0^1]- (This 
includes the case of the FPU problem.) The slow energy exchange behavior of this system was 
analyzed in detail by |Hairer et al.|(2006} XIII, see especially p. 495) using the so-called modulated 



Fourier expansion; we now give a brief, high-level summary of this work. In the notation of Hairer 
et al. the exact solution is asymptotically expanded as x*(f) ~ Cfit), where 

x,[t) = yit] + e''''zit) + e-''''z[t). 



*Note that although Step 2 might appear to be ill-defined, due to the fact that fl is possibly singular, the singularity can 
be removed by substituting the relation hD./2 = tan (^hD./2^ . The explicit equation for q„+i and is then calculated 
to be 

qn+A_f cos(hn) h/2[l + cos{hn))\ f qA 

Pn+iJ ~ [-2/h (1 - cos [hn]) cos [hn) j[ptj' 

which is seen to recover the correct, purely kinematic equation when n = Q = 0. 
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Here, y (?) is a smooth real-valued function and z[t) is a smooth complex-valued function, and 
these can be partitioned according to the blocks of f2 as y = (jo-yi) and z = (zq.Zi). They 
show that, assuming the exact solution has energy bounded independent of co, this implies 
z{t) = [oj~^), so z describes the slow-scale evolution of the system. Plugging in this ansatz for a 
highly oscillatory system, and eliminating the variables yi and zq, the slow evolution turns out to 
be described by 

Zicozi = (yo,0) zi + O (oj-^) , 

where here the slow force g[q) = —VU[q) has also been block-decomposed as g = (go> gi]- 

Hairer et al.| compare the above with the numerical solution for a trigonometric integrator, 
which is similarly expanded as 

Xh{t) = yh{t) + e''''zh{t) + e~''''zh{tl 

with yii = [yh,o>yh,i] and z;, = [zh,o,Zh,i]- For the unmodified Deuflhard/impulse method, in 
particular, the slow-scale numerical evolution is given by 

which implies that the equation for the impulse solution is consistent with that for Z\. 

We now finally have what we need to prove our main result on the slow energy exchange 
behavior of the variational IMEX method. 



Theorem 4.1. Let the variational IMEX method he applied to the highly oscillatory problem above, 
and suppose the numerical solution remains bounded. Then the ordinary differential equation for 
Zh,i, describing the slow energy exchange in the numerical solution, is consistent with thatforzi in 
the exact solution; this holds up to order O {co~^] ■ 



Proof As we have previously shown, the IMEX scheme corresponds to a modified impulse method 
with frequency o). Therefore, to get the equation for we must modify the equation above 
for z/j i by replacing o) with o) on the left hand side. However, notice that Step 2 of the modified 
method advances the original state vector j , rather than the modified ^ j . Changing from 

[^^" ) to {j'-^^ j also introduces a scaling factor of ca/oj on the right hand side. Therefore, the 
variational IMEX solution satisfies the slow-scale equation 

2ia)2h,i = — -I — [yh.o.O] Zh,i + O {co~'^(b] . 

CO OXi ^ ' 

Finally, cancelling the (b factors and multiplying by cu, we once again get 

2icolh,\ = iyh,o,0) Zh,i + O (w"^) , 

which is the same as that for the coefficient z^^i in the original, unmodified impulse method. This 
completes the proof. □ 



16 



A. STERN AND E. GRINSPUN 



Acknowledgments 

The authors would like to thank Will Fong and Adrian Lew for helpful conversations about 
the stability of AVI methods, Houman Owhadi for suggesting that we examine the FPU problem, 
Dion O'Neale for advice regarding the FPU simulation plots, Teng Zhang for pointing out a bug 
in an earlier version of our FPU simulation code, and Jerry Marsden and Reinout Quispel for 
their valuable suggestions and feedback about this work. We especially wish to acknowledge 
Marlis Hochbruck, Arieh Iserles, and Christian Lubich for suggesting that we try to understand 
the variational IMEX scheme as a modified impulse method; this perspective led us directly to 



the results in Section 4 Finally, this paper benefited greatly from the thoughtful critiques and 



suggestions of the anonymous referees, whom we also wish to thank. 

References 

Ascher, U. M., and S. Reich (1999a), The midpoint scheme and variants for Hamiltonian systems: 
advantages and pitfalls. SIAMJ. Sci. Comput., 21 (3), 1045-1065 (electronic). 

Ascher, U. M., and S. Reich (1999b), On some difficulties in integrating highly oscillatory Hamilton- 
ian systems. In Computational molecular dynamics: challenges, methods, ideas (Berlin, 1997), 
volume 4 of Lecture Notes in Computational Science and Engineering, pages 281-296. Springer, 
Berlin. 

Ascher, U. M., S. J. Ruuth, and R. J. Spiteri (1997), Implicit- explicit Runge-Kutta methods for 
time-dependent partial differential equations. Appl. Numer. Math., 25 (2-3), 151-167. Special 
issue on time integration. 

Biesiadecki, J. J., and R. D. Skeel (1993), Dangers of multiple time step methods. /. Comput. Phys., 
109 (2), 318-328. 

Boxerman, E., and U. Ascher (2004), Decomposing cloth. In SCA '04: Proceedings of the 2004 ACM 
SIGGRAPH/Eurographics Symposium on Computer Animation, pages 153-161. Eurographics 
Association, Aire-la-Ville, Switzerland, doi : 10 . 1145/1028523 . 1028543, 

CrouzeiK, M. (1980), Une methode multipas implicite-explicite pour 1' approximation des equa- 
tions d'evolution paraboliques. Numer Math., 35 (3), 257-276. 

Deuflhard, P. (1979), A study of extrapolation methods based on multistep schemes without 
parasitic solutions. Z. Angew. Math. Phys., 30 (2), 177-189. 

Eberhardt, B., 0. EtzmuE, and M. Hauth (2000), Implicit- explicit schemes for fast animation with 
particle systems. In Proceedings of the 11th Eurographics Workshop on Computer Animation 
and Simulation (EGCAS), pages 137-151. Springer, Interlaken, Switzerland. 

Farr, W. M., and E. Bertschinger (2007), Variational integrators for the gravitational A/^-body 
problem. Astrophys. J., 663 (2), 1420-1433. arXiv : astro-ph/0611416, 

Fermi, E., J. Pasta, and S. Ulam (1955), Studies of nonlinear problems. Report LA-1940. Los Alamos 
National Laboratory, Los Alamos, NM. 

Fong, W, E. Darve, and A. Lew (2007), Stability of asynchronous variational integrators. 
In PADS '07: Proceedings of the 21st International Workshop on Principles of Advanced 
and Distributed Simulation, pages 38-44. IEEE Computer Society Press, Washington, DC. 
[doi:10.110 9/PADS .2007.29, 



IMPLICIT- EXPLICIT VAHIATIONAL INTEGRATION OF HIGHLY OSCILLATORY PROBLEMS 



17 



Galgani, L., A. Giorgilli, A. Martinoli, and S. Vanzini (1992), On the problem of energy equipartition 
for large systems of the Fermi-Pasta-Ulam type: analytical and numerical estimates. Phys. D, 
59 (4), 334-348. 

Hairer, E., C. Lubich, and G. Wanner (2003), Geometric numerical integration illustrated by the 
Stormer-Verlet method. Acta Numer., 12, 399-450. 

Hairer, E., C. Lubich, and G. Wanner (2006), Geometric numerical integration, volume 31 
of Springer Series in Computational Mathematics. Springer- Verlag, Berlin, second edition. 
Structure-preserving algorithms for ordinary differential equations. 

Izaguirre, J. A., Q. Ma, T. Matthey J. Willcock, T. Slabach, B. Moore, and G. Viamontes (2002), Over- 
coming instabilities in Verlet-I/ r-RESPA with the mollified impulse method. In T. Schlick and 
H. H. Gan, editors, Computational Methods for Macromolecules: Challenges and Applications — 
Proceedings of the 3rd International Workshop on Algorithms for Macromolecular Modeling, 
volume 24 of Lecture Notes in Computational Science and Engineering (LNCSE), pages 146-174. 
Springer, Berlin. 

Leimkuhler, B. J., S. Reich, and R. D. Skeel (1996), Integration methods for molecular dynamics. In 
Mathematical approaches to biomolecular structure and dynamics (Minneapolis, MN, 1994), 
volume 82 of IMA Vol. Math. Appl, pages 161-185. Springer, New York. 

Lew, A., J. E. Marsden, M. Ortiz, and M. West (2003), Asynchronous variational integrators. Arch. 
Ration. Mech. Anal, 167 (2), 85-146. 

Marsden, J. E., and M. West (2001), Discrete mechanics and variational integrators. Acta Numer, 
10, 357-514. 

McLachlan, R. I., and D. R. J. O'Neale (2007), Comparison of integrators for the Fermi-Pasta-Ulam 

problem. Preprint NI07052-HOP. Isaac Newton Institute for Mathematical Sciences, Cambridge, 

UK. Available from: |http: //www.newton. ac .uk/preprints/NI07052 .pdf ( 
McLachlan, R. I., and G. R. W. Quispel (2002), Splitting methods. Acta Numer, 11, 341-434. 
Moser, J., and A. P. Veselov (1991), Discrete versions of some classical integrable systems and 

factorization of matrix polynomials. Comm. Math. Phys., 139 (2), 217-243. 
Suris, Y. B. (1990), Hamiltonian methods of Runge-Kutta type and their variational interpretation. 

Mat. Model., 2 (4), 78-87. 
Zhang, G., and T. Schlick (1993), LIN: a new algorithm to simulate the dynamics of biomolecules 

by combining implicit-integration and normal mode techniques. /. Comput. Chem., 14 (10), 

1212-1233. doi: 10. 1002/jcc.540141011, 
Zhang, M., and R. D. Skeel (1997), Cheap implicit symplectic integrators. Appl. Numer Math., 25 

(2-3), 297-302. Special issue on time integration. 

Ari Stern, Department of Applied and Computational Mathematics, California Institute of Technology, 
Pasadena, California 91125 

Current address: Department of Mathematics, University of California, San Diego, La JoUa, California 92093-01 12 
E-mail address: astern@math.ucsd. edu 



Eitan Grinspun, Department of Computer Science, Columbia University, New York, New York 10027 
E-mail address: eitan@cs . Columbia, edu 



